Pairing mean-field theory for the dynamics of dissociation of molecular Bose-Einstein condensates 
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We develop a pairing mean-field theory to describe the quantum dynamics of the dissociation of molecu- 
lar Bose-Einstein condensates into their constituent bosonic or fermionic atoms. We apply the theory to one, 
two, and three-dimensional geometries and analyze the role of dimensionality on the atom production rate as a 
function of the dissociation energy. As well as determining the populations and coherences of the atoms, we cal- 
culate the correlations that exist between atoms of opposite momenta, including the column density correlations 
in 3D systems. We compare the results with those of the undepleted molecular field approximation and argue 
that the latter is most reliable in fermionic systems and in lower dimensions. In the bosonic case we compare 
the pairing mean-field results with exact calculations using the positive-P stochastic method and estimate the 
range of validity of the pairing mean-field theory. Comparisons with similar first-principle simulations in the 
fermionic case are currently not available, however, we argue that the range of validity of the present approach 
should be broader for fermions than for bosons in the regime where Pauli blocking prevents complete depletion 
of the molecular condensate. 

PACS numbers: 03.75.-b, 03.65.-w, 05.30.-d, 33.80.Gj 



I. INTRODUCTION 

The generation and detection of strongly correlated atomic 
ensembles is becoming one of the central themes in the study 
of ultracold quantum gases, leading to the birth of the subfield 
of quantum-atom optics. Recently there have been an increas- 
ing number of experiments probing the higher-order coher- 
ences and correlation functions that are fundamental to this 
field and are analogous to similar developments in the early 
days of quantum optics with photons. 

One of the first experiments in this area was the measure- 
ment of a two-time second-order correlation function in a 
thermal gas of metastable Neon atoms falling on four mi- 
crochannel plate (MCP) detectors UJ. Another early exper- 
iment was the measurement of the rate of three-body losses 
from a Bose-Einstein condensate (BEC), from which the lo- 
cal third-order correlation function could be inferred (see 
also iU). More recently, the arsenal of experimental tech- 
niques to study higher-order correlations has been expanded to 
include the measurement of photoassociation rates M, shot- 
noise spectroscopy of absorption images IH |^ IH |8[], atom 
counting using fluorescence imaging and high-finesse optical 
cavities ll^O, and position-resolved counting of metastable 
Helium atoms using arrays of MCP detectors [ 111 Il2l iTsll . 
Proposals for additional techniques include atom counting us- 
ing photoionization [1 14!], and freeze-in techniques with optical 
lattices for in-situ measurements of non-local spatial correla- 
tions L15i1 . 

The experimental systems that are being studied with these 
new techniques for correlation measurements include Bose- 
Einstein condensates and degenerate Fermi gases ii[ni[ii, 

atom lasers lEflL one-dimensional (ID) and quasi- ID Bose 
gases JH 3, Bose and Fermi gases in optical lattices 

la 01, four- wave mixing via BEC collisions |13] and molec- 
ular dissociation of BECs |j5|l. The last two systems are rel- 
evant to the present work and have similar atom-atom corre- 



lations on the dissociation sphere and on the s-wave scatter- 
ing halo, respectively 1 16, 17J. All the experiments described 
above generate quantum correlations that cannot be described 
by Gross-Pitaevskii mean-field theory. They provide a chal- 
lenge for theorists to develop and apply more advanced tech- 
niques in order to quantitatively model experimental systems, 
such as used in Refs. Ill [11 [11 113, III 111 111 Sill 111]. 

Particular interest in molecular dissociation experiments 
stems from the prospect of gen erating and dete cting squeezed 
states of matter waves ITslTTl IH |M IH, HI HI, with po- 
tential applications in precision measurements. There also 
exists the possibility of producing Einstein-Podolsky-Rosen 
correlations and other entangled states of massive particles 
|32, 33, 34, 35, 36, 37, 38, 39] for fundamental tests of quan- 
tum theory in mesoscopic regimes and for quantum informa- 
tion science applications. Additionally, molecular dissocia- 
tion can serve as a probe of two-body interactions, collisional 
resonances, and spectroscopic properties of Feshbach reso- 
nance molecules 1401 l4ll l42l |43|] . Other recent studies of 
dissociation are concerned with the role of confinement on 
the stability of a molecular BEC against dissociation |44], the 
effect of magnetic field fluctuations on dissociation |45,], and 
the dynamics of molecule-atom conversion in optical lattices 



The purpose of this paper is to develop and apply a pairing 
mean-field theory (PMFT) to the dynamics of dissociation of a 
BEC of molecular dimers consisting of pairs of either bosonic 
or fermionic atoms. The mean fields for the molecules are 
introduced in the standard way at the level of amplitudes, 
whereas for the atoms they are introduced at the level of pair- 
ing fields Ii47i1 that describe normal and anomalous moments 
of the atomic creation and annihilation operators. This ap- 
proach takes into account the molecular depletion and can fa- 
cilitate the calculation of atom-atom pair correlation functions 
beyond the coherent level of the simple mean field theory. The 
present PMFT is similar to that of Jack and Pu Ii23il . except 
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that: (/) we extend the analysis to one-dimensional (ID) and 
two-dimensional (2D) systems; (//) we relax the approxima- 
tion of the constant density of atomic states around the peak 
dissociation energy; and (///) we calculate the pair correlation 
functions for dissociated atoms. 

First-principle quantum simulations of the dissociation of 
molecular BECs into bosonic atoms have been performed us- 
ing the stochastic positive-P representation method |18, 
I24I1 . however this technique has some limitations. Most signif- 
icant is the short time scale over which the positive-P method 
remains valid and produces results free of large sampling er- 
rors or the boundary term problem |48]. For pure dissociation, 
without s-wave scattering interactions, the typical timescale 
for the positive-P simulations is limited to about 50% conver- 
sion 1 24]. The addition of s-wave scattering interactions fur- 
ther limits this to even shorter durations and only 5% — 10% 
conversion for typical experimental parameters ll24ll . The sit- 
uation is less clear in the fermionic case; even though similar 
first-principle techniques for simulating multi-mode fermion 
dynamics have recently been developed |49], these techniques 
have not been applied yet to the problem of multi-mode dis- 
sociation for quantitative comparison. 

Given the limitations of exact techniques, it is useful to de- 
velop alternative, approximate theoretical techniques that can 
deal with longer durations of dissociation and larger conver- 
sion percentage. The PMFT developed here and in Ref. ll23ll 
falls into this category and goes beyond previous approximate 
theories that ignore molecular depletion altogether ll24l Isoll . 
As such, the range of its validity with respect to molecular 
depletion covers dissociation durations corresponding to near 
complete conversion. In addition, the PMFT is able to treat 
the dissociation into bosonic and fermionic atoms on equal 
footing and therefore is closer to being able to quantitatively 
describe the experiments on dissociation of ^°K2 dimers into 
fermionic atoms 1|5|], in addition to the earlier experiments on 
dissociation of ^'^Na2 ifsil l and ^^Rb2 ll52ll dimers into bosonic 
atoms. 

This paper is organized as follows. In Sec. |II]we introduce 
the model Hamiltonians describing the dissociation into dis- 
tinguishable and indistinguishable atoms and discuss the va- 
lidity of approximations involved in our model. In Sec. Ulllwe 
give the simplest possible description of dissociation based on 
a Fermi's golden rule calculation of the molecular decay rate 
and discuss the implications for dissociation in ID, 2D and 
3D geometries. In Sec. |IV]we develop the PMFT and discuss 
our main results for dissociation dynamics, atomic momentum 
distribution, and atom-atom correlations. Our conclusions are 
summarized in Sec. |VT] 



II. THE MODEL 

In our treatment of dissociation we assume an initial con- 
dition of a stationary BEC of diatomic molecules consisting 
of either (/) two distinguishable fermionic or bosonic atoms 
(e.g. in two different spin states) or (//) two indistinguishable 
bosonic atoms. In the case (/), we further assume that both 
constituents are of the same species and the same mass; we do 



not consider heteronuclear molecules. The model Hamiltoni- 
ans for each of these situations are described below, and we 
discuss the initial conditions and the validity of our approxi- 
mations in Sec. Ill CI 

A. Dissociation into distinguishable atoms 

The quantum field theoretic effective Hamiltonian describ- 
ing the system in one, two, or three (D — 1,2,3) spatial di- 
mensions is given, in a rotating frame, by Ii53i1 

-ihxD (*J*1*2 - ^^^-I^o) } ■ (1) 

The molecular field is described by a bosonic operator 
\['o(x, i) whilst the two atomic fields are described by ei- 
ther bosonic or fermionic operators ^i(x, i) and ^2(^1 
The field operators satisfy the respective commutation or 
anti-commutation relations, i.e., [^'^(x, <), >I'j(x', t)] = 
S^JS'•^\x - x') in the case of bosonic atoms or 
{^',(x,i),^'](x',i)} = %5(^)(x-x') {I = 1,2) in the 
fermionic case. For notational simplicity, we use x for the 
position in ID, 2D and 3D cases, with the understanding that 
in the ID case it is a scalar. 

The first term in the Hamiltonian (|T) describes the kinetic 
energy where the atomic masses are mi = 7712, whereas the 
molecular mass is toq = mi +7712 ~ 2mi. The coupling term 
Xd is responsible for coherent conversion of molecules into 
atom pairs, e.g. via optical Raman transitions or a Feshbach 
resonance sweep (see, for example, Refs. ll47[ Issl l54. 55, 5^ 
UlllSSt] and ll59il60ll for recent reviews). For the Raman case, 
the coupling xzd is expressed in terms of the Rabi frequencies 
for free-bound and bound-bound transitions as in Ref. llssll . 

In the case of a Feshbach resonance, the coupling XciD is 
given by (see [gTlllilll for notational consistency) 



X3D ^ y 4:TratgA^AB/mi. (2) 

Here Afi — i^i + IJ.2 — Mo is the difference in the magnetic 
moments of the atomic and the bound molecular states, AP 
is the magnetic width of the resonance, and abg is the back- 
ground scattering length for s-wave collisions of the atoms in 
the two spin states. In systems of reduced dimensionality (ID 
or 2D) and away from confinement induced resonances, the 
couplings xiD and X2D, are obtained by integrating over the 
ground-state wave function in the tightly confined dimensions. 
Assuming harmonic trapping potentials in the eliminated di- 
mensions, with oscillation frequencies lo±_ = ujy = uj^ in the 
ID case and lOz in the 2D case, the respective atom-molecule 
couplings are 

XiD = X3d/ (3) 
X2D = xw/{2nllf'\ (4) 
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Here l± = y^h/miujj_ and Iz = ^JTiJm^Lol are the harmonic 
oscillator lengths for the atoms. For a given transverse trap 
frequency wj, in the ID geometry, the dissociation detuning 
must satisfy |A| < o^x to avoid transverse excitations and 
maintain the validity of the ID treatment. Similarly, in the 
2D case one has to satisfy | A| < in order that the system 
remains in the 2D regime in the xy plane. 

The quantity A in the Hamiltonian ([Hi is the effective de- 
tuning between the molecular state with energy Eq and the 
atomic states with energies and E^- It is defined differ- 
ently depending on whether the coupling is due to a sweep 
through a Feshbach resonance, or due to Raman photoassoci- 
ation lasers. 

(0 In the Feshbach resonance case, dissociation of an ini- 
tially stable (So < i?i + S2) molecular BEC into the 
constituent atoms may be achieved by a rapid mag- 
netic field sweep onto the atomic side of the resonance 
(i?o > E\+ Ei) [52]. The detuning is defined as the 
overall energy mismatch 2ftA = Ei + E2 ~ Eq be- 
tween the free two-atom state in the dissociation thresh- 
old with energy Ei + E2 and the bound molecular 
state of energy Eq. The sweep is implemented to re- 
sult in a finite and negative detuning 2hA < after the 
sweep {Eq > Ei + E2), in which case the molecules 
become unstable against spontaneous dissociation into 
free atom pairs. 

(ii) In the case of two-photon Raman photoassociation fl^, 
the overall energy mismatch is given by 2hA ~ 
El + E2 — Eq — huj, where u is the frequency differ- 
ence between the two Raman lasers. Similarly, for an rf 
transition |5] is the rf frequency. 

The Hamiltonian ([T]l conserves the total number of atomic 
particles 

N = 2No{t) + Ni{t) + N2{t) ^ const, (5) 

where the constant is given by 27Vo(0) for the vacuum initial 
condition for the atoms. 



B. Dissociation into indistinguisliable atoms 

Molecular dissociation into pairs of indistinguishable 
bosonic atoms in the same internal state is described by the 
following Hamiltonian ifsill . in a rotating frame 




_,Y(*s*?-*r*o)}. (6) 

where ^1 (x) is the atomic field operator, and xd is the atom- 
molecule coupling. The microscopic expression for the 3D 



coupling X3D in the case of a Feshbach resonance is given by 



X3Z5 = ySTrflbgA^AB/mi, (7) 

where abg is the atom-atom background scattering length, 
Afi = 2fii — /io is the magnetic moment difference be- 
tween the atomic and the bound molecular states, and AB 
is the magnetic width of the resonance. For systems of re- 
duced dimensionality, Eqs. Q and (|4]i for xiD and X2D are 
unchanged. 

The Hamiltonian (|6]l conserves the total number of atomic 
particles 

N = 2No{t) + Ni{t) = const, (8) 

where the constant is given by 27Vo(0) for the vacuum initial 
condition for the atoms. 



C. Initial conditions and approximations 

Starting with a stable molecular condensate, we assume that 
the coupling xd is switched on in the regime of a sudden jump 
ll43ll : initially the molecules are assumed to be in a coherent 
state, whereas the atomic fields are in the vacuum state. The 
energy level configuration after the Feshbach sweep or the Ra- 
man transitions (corresponding to A < 0) is the actual initial 
condition for our simulations. 

For a molecule at rest, the excess of potential energy 2^| A| 
(which we also refer to as the dissociation energy) is converted 
into kinetic energy 2?i^fc^/ (2mi) of atom pairs with equal but 
opposite momenta ±ko, where fcp — |ko| — y/2mi\A\/h. 
This is the physical origin of the expected correlations be- 
tween the atoms, which we will study below in the context 
of many-body field theory. 

The main limitation of our treatment is that we consider a 
spatially uniform system in a cubic box with periodic bound- 
ary conditions. The question that one has to address then is 
how well the results of a uniform model can describe realis- 
tic nonuniform systems. It has previously been shown ll24ll 
that while the uniform treatment can give quantitatively rea- 
sonable results for the total number of atoms and the atomic 
density distribution (provided the uniform system is appropri- 
ately size-matched to a nonuniform system, as in Sec. |IV All ), 
such a treatment is not adequate for giving correct quantitative 
results for density-density correlation functions. The reason 
is the mode-mixing due to inhomogeneity of trapped conden- 
sates, which can strongly degrade the correlations compared 
to the predictions of the uniform model |24]. Nevertheless, 
the merit of calculating the correlation functions here is that 
the results provide upper bounds for the correlations |24, 5C|] 
and hence aid qualitative understanding. 

Another obvious consideration when applying the uniform 
results to size-matched nonuniform systems is that the dura- 
tion of dissociation should not exceed the time required for 
the atoms to propagate distances larger than the size of the 
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uniform box L. Thus the resuhs are valid for times 



t<t 



L 



max - — - iV"?l/(2ft|A|), 

Wo 



(9) 



where vq — hko/nii is the mean velocity of the dissociated 
atoms. In other words, while our numerical results can be for- 
mally obtained for arbitrarily long times, the above equation 
should be used to cutoff these results at tmax when they are 
applied to a specific nonuniform system of size L. 

The next major assumption in our treatment is the mean- 
field coherent state for molecules at all times lEsll . This ap- 
proximation breaks down once the dissociation approaches 
the regime of complete depletion of the molecular condensate 
(~ 100% conversion) when quantum fluctuations become in- 
creasingly important. Accordingly, the results of the mean 
field theory can not be trusted past that point in time; we 
present these results only for academic purposes. We note that 
nearly complete conversion in the mean field theory would al- 
ways take place in the case of bosonic atoms, but not neces- 
sarily in the fermionic case where the Pauli exclusion princi- 
ple for the atoms can dominate the molecular depletion. In 
this situation the fermionic results can prove to be reliable for 
longer durations than in the bosonic case. 

Finally, we do not include any s-wave scattering interaction 
terms. Potentially these interactions can affect the atom-atom 
correlations in realistic nonuniform systems, especially in the 
long time limit. However, their effect is negligible at low par- 
ticle densities and short dissociation times |24]. At the level 
of mean fields, the s-wave scattering terms can be neglected 
if the dissociation energy is much larger than the mean-field 
interaction energy. More importantly, it has been shown in 
Ref. [24] that the strongest degradation of atom-atom cor- 
relations in realistic systems (compared to the predictions of 
uniform models) comes not from the s-wave scattering in- 
teractions, but from the mode-mixing due to inhomogeneity. 
Accordingly, any theory that attempts to give quantitatively 
correct results for atom-atom correlations has to first address 
the question of inhomogeneity before including s-wave inter- 
actions. The pairing mean field theory developed here does 
not have this goal, but only represents an intermediate step 
between the simple analytic theory of undepleted molecules 
ll24l [soil and a more complete theory that can reliably treat 
mode-mixing due to inhomogeneity and s-wave scattering in- 
teractions for both bosonic and fermionic atoms. 



m. DECAY RATE AND DISSOCIATION DYNAMICS 
FROM FERMI'S GOLDEN RULE 



regime. In particular, Fermi's golden rule gives the correct 
initial behavior of the total atom number as a function of time 
since at this stage neither bosonic stimulation nor Pauli block- 
ing affects the dynamics of individual modes. It also provides 
simple and useful insights into the differences associated with 
the dimensionality of the system, which are not explicit in 
the numerical results. The derivation presented here is a di- 
rect calculation based on the actual many-body Hamiltonian 
rather than an estimate based on the comparison of the mean- 
field energy shift and the energy shift of a Feshbach resonance 
state as in Ref. Isill . 

According to Fermi's golden rule, the molecular decay rate 
r is given by L5Li61J 



2-K 

r=—\V,na\''D 



(2) 



(10) 



where £'(^'(e) is the density of two-atom states, e is the total 
kinetic energy of two atoms corresponding to the total disso- 
ciation energy e = 2/i|A|, and Vma is the transition matrix 
element for the atom-molecule coupling channel. Calculat- 
ing Vma explicitly from the interaction Hamiltonian in Eq. ([T]) 
(see Appendix IaI for details), we obtain the following results 
in 1,2, and 3 dimensions: 



r = < 



^ X2Z5' 



(ID), 
(2D), 



(11) 



[ 2A3xL|ArV^, (3D) 



A. Distinguishable atoms 



where A = (mi/2fi,)^/^ is a constant. Since each molecule 
produces two atoms (one per spin state) we can treat the decay 
rate T as the atom production rate using the conserved total 
pai-ticle number 2iVo(t) + Ni{t) + N2{t) = 27Vo(0) = const. 
Here, N.i = {Ni) with -/Vo(O) being the total initial number of 
molecules with no atoms present. Thus, the solution iVo [t) = 
Afo(O) exp(— rt) to the rate equation for the molecular decay, 
can be rewritten as iVi(t) + Ar2(t) = 2iVo(0) [1 - exp(-rt)]. 
Noting also that iVi {t) — N2 (t) since iVi — N2 is conserved, 
the final result for the total number of atoms in each spin state 
takes the following simple form: 

N^2){t) = Nomi - e-^') ^ No{0)Tt, (H « 1). (12) 

We see immediately that the atom production rate in ID, 2D 
and 3D increases, remains unchanged and decreases, respec- 
tively, with the dissociation energy 2fi,|A|, following the re- 
spective dependence of the decay rate F on the detuning | A|. 
This is a direct consequence of the different dependence of the 
density of atomic states D{E) on energy E in ID, 2D and 3D. 



Before presenting our numerical results within the frame- 
work of the PMFT, we discuss the simpler results that follow 
from a Fermi's golden rule calculation of the molecular decay 
rate. Due to the conservation of the total number of atomic 
particles, the molecular decay rate can be converted into the 
rate of atom production and serve as the simplest descrip- 
tion of the dynamics of dissociation in the initial spontaneous 



B. Indistinguistiable atoms 

The decay rate F derived above is for dissociation into dis- 
tinguishable atoms. For completeness, we now present the 
derivation in the case of indistinguishable (bosonic) atoms and 
show that the resulting expression in 3D coincides with the ex- 
pression given in Ref. |51]. The reasons for presenting these 
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details are two-fold. First, we would like to point out that even 
though the present paper is primarily concerned with dissocia- 
tion into distinguishable atom pairs (either fermions or bosons 
in different spin states), the results in the bosonic case can be 
easily modified to describe the case of indistinguishable atoms 
in the same spin state. Secondly, we would like to spell out the 
connection of the microscopic Feshbach parameters with the 
coupling constant X3D in the respective field-theory Hamil- 
tonians without leaving room for confusion about "factors of 
two". The details of derivation of the decay rate F' in the 
present case are given in AppendixiBl the final result is: 



The Hamiltonian of the system, Eq. ([T]i, can now be written 



r Ax?,,/2|A|i/2, (ID), 



(2D), 



(13) 



_ A^xLlAr/V'T, (3D) 



For the same coupling strength xd these rates appear to be 
simply half the corresponding result for the distinguishable 
case, Eq. (fTTl i. However, one has to take into account that the 
interaction Hamiltonian for the present indistinguishable case 
has an extra factor of 1/2 in front of xd and that the definition 
of Xd in the two cases are different. 

The total number of atoms produced during the initial stage 
of dissociation (T't <C 1) is 



iVi(i) = 2iVo(0)(l- 



-r't\ 



2No{0)T't, (F'i < 1). (14) 



Comparing this with the result of Eq. ( fT2] l we see that for the 
same coupling constant xd in the Hamiltonians (|6]l and ([T]), 
the total number of atoms produced in the indistinguishable 
case corresponds to the number of atoms in each spin state 
produced in the distinguishable case. 



IV. PAIRING MEAN-FIELD THEORY 

A. Distinguishable atoms 

Considering a uniform system in a cubic box of volume 
V ^ (or a square of area A = in 2D, or a line of length 
L in ID), we expand the molecular and atomic field operators 
in terms of plane-wave modes 

*o(x,t) = _1^^^6k(^)e^''■^ (15) 

^j(^'*) = IS72Ek°J^k(i)e*''", (16) 

where 6k (i) and aj_k(i) ij = 1,2) are the corresponding 
creation operators satisfying the usual bosonic commutation 
relation for molecules [Sk, S]^,] = (5k. k' and fermionic anti- 
commutation or bosonic commutation relations for the atoms. 

Here 



ai,k, a] ,^.,} = (5k,k"5i,i or [0,1^,0,] ^^,] = (5k,k"5i 



k 
n 



-'1,3- 



(27r/L)n is the momentum in wave-number units, with 
riz) and rii = 0, ±1, ±2, . . .(i = x,y, z) in the 



as 



E 

k,i=l,2 



2m,- 



E(^l"i.ka2,q- k - al q_^a\ J>q). (17) 



n. 

To proceed, we next assume that the initial condition is a 
large molecular condensate in the mode 60 ((&o(0)^o(0)) ^ 
1) and that the dynamics of dissociation retains the molecular 
population in the same condensate mode 60 (0- This implies 
that we ignore any atom-atom recombination process in which 
two atoms with non-opposite momenta ai k and a2.q-k com- 
bine to populate a non-condensate molecular mode 6q with 
q 7^ 0. This is a reasonable approximation at least during the 
initial stages of dissociation when recombination is negligible. 
Thus, the Hamiltonian ( [TtI i now reduces to 

H = y^^_ ^Ak fal, kQi,k + fl2.k°2,k) 

-iftK^^(6jai^ka2 -k - a^-k^lk^o), (18) 

where we have introduced Ak = A + /ik^/ (2toi) and a new 
coupling constant 



XD 



(19) 



3D case, with similar relations in the ID and 2D cases. 



We note that the Hamiltonian ( fTsT i conserves the total number 
of atomic particles, N — 26q6o + J2k ("-i k + ?^2.k), and the 
number difference, rii k — "■2,k, where n^-^k = a]k'^j k (j — 
1, 2) are the particle number operators for the atoms. 

The next step in the mean-field approach is to assume that 
the molecular condensate is initially in a coherent state with 
an amplitude /3o ((6^,(0)60(0)) = |/3oP > 1), and that 
it remains in a coherent state during the dynamical evolu- 
tion. The condensate dynamics are then treated at the level 
of the mean-field amplitude, /3{t) — {bo{t)), which implies 
that we can approximate the higher-order correlation function 
(6ja1.ka2.-k) ~ /3o(ai,ka2,-k)- Thus, the decorrelation as- 
sumption is imposed at all times, which is the main limitation 
of the method. This treatment will become increasingly inad- 
equate as the molecular population depletes and approaches 
zero. 

Once we impose the decorrelation assumption, the dynam- 
ics of the atomic fields can be described in terms of the normal 
and anomalous populations 

nu{t) = {ni,k{t)) = (n2,k(i)), (20) 
rukit) = (mk(t)) = (ai,k(i)a2,-k(i)), (21) 

which describe the occupation numbers of the atomic modes 
with momenta k, Eq. ( |20l i, and the correlation between the 
modes with equal but opposite momenta, k and — k, in the 
two spin states, Eq. ( l2Tl i. In this description, higher-order cor- 
relation functions factorize according to Wick's theorem, i.e., 
they can be expressed in terms of the second-order moments 
nk(i) and m^{t). 
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Writing down the Heisenberg operator equations of motion 
following from the Hamiltonian ( fTSl l and imposing the pairing 
mean-field decorrelation approximation, one can show that 
the equations of motion for the scaled mean-field amplitude, 



fir) = p{t)/f3o, 



(22) 



and the normal and anomalous populations n\^{t) and m^{t) 
form a closed set and can be written as 



dr 
dr 

dfjr) 
dr 



= mliT)f{T)+m^{T)r{T), 

= -2iJk™k(r)+/(T)[l±2nk(T)], (23) 
1 , 



mk(T). 



In these equations and throughout the rest of this paper, the 
upper (lower) sign refers to bosonic (fermionic) atom statis- 
tics, and we have introduced dimensionless time t = t/to and 
dimensionless effective detuning 



- Akto - [?ikV(2mi) + A]to = 



(24) 



where to — 1/(k/3o) is the time scale (where we have as- 
sumed that /3o is real without loss of generality), g = |q| is 
the absol ute value of the dimensionless momentum q = kdo, 
do = ^JhtQ/{2mi) is the length scale, and 5 ~ Aip is the 
dimensionless bare detuning. The quantity /(r) is the frac- 
tional molecular- ampHtude so that |/(t)P — Nq{t)/Nq{Q) 
corresponds to the fraction of molecules relative to their 
initial total number Equations ( |23] | conserve the quantity 
|'Tik(^)P — f^k(i)[l ± n-i^it)] as its time derivative is zero, 
and so the pair and anomalous densities are directly related 
according to 



|mk(i)P =nk(i)[l±«kW]. 



(25) 



This will be useful in the calculation of correlation functions 
in Sec.rVB] 

Equations ( |23T l are equivalent to those solved in Ref. ||23ll 
for describing dissociation dynamics in 3D using spin- 1/2 
Pauli matrices |62]. Here, we extend this treatment to ID 
and 2D systems. In addition, we analyze the second-order 
correlation functions for the atoms and give detailed quanti- 
tative assessment of the approximations involved in the pair- 
ing mean-field method compared to the exact first-principles 
treatment using the positive P-representation for bosons i24ll . 



1. Results for total atom numbers 

In Figs. [T] and |2] we show typical results from solving 
Eqs. ( |23] ) for fermionic and bosonic atom statistics respec- 
tively. We plot the fraction of the total number of atoms in 
each spin component relative to the total initial number of 
molecules, A^i,2(''')/-/Vo(0), as a function of time. The two 
graphs, (a) and (b), are for two different values of the dimen- 
sionless detuning 5 — Aip {5 < 0) and the same initial dimen- 
sionless molecular density NoiO)/l^ = 3.1 (D = 1,2,3), 
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FIG. 1 : Dissociation of a molecular condensate into fermionic atoms. 
We plot the total number of atoms in each spin component A^i,2(t) 
[Ni{t) = A^2('r)] relative to the total initial number of molecules 
A'^o(O) as a function of the dimensionless time r — t/to, where to = 
l/hi\/ No{0) is the time scale. The graphs (a) and (b) are for |<5| = 16 
and \ 5\ = 3174, respectively. The different curves are for ID, 2D and 
3D cases as shown. In all cases, the initial dimensionless molecular 
density is No{0)/1° = 3.1 (D = 1, 2, 3), where I = L/do is the 
dimensionless length. The dotted lines next to each solid line are the 
results from Fermi's golden rule, Eqs. ( jilt and il2\ . See text for 
further details. 



where I — L/do is the dimensionless length corresponding to 
the quantization length L, with dg — ^yhto/{2mi) being the 
length scale. 

These parameters can be size-matched with a realistic 
nonuniform system corresponding, for example, to a BEC 
of molecular dimers made of fermionic "^^K atoms as in 
Ref. [5]. Our empirically derived prescription for the size- 
matching 1240 is as follows. Assuming for simplicity an 
isotropic harmonically trapped molecular condensate with a 
Thomas-Fermi parabolic density profile, we match the pa- 
rameters of the present uniform treatment to have the same 
peak density po as the trapped condensate and choose the 
length L of the uniform box to result in the same ini- 
tial total number of molecules as in the trapped conden- 
sate Nq{0) = (87r/15)i?j.^po = L^Po, where Rtf is the 
Thomas-Fermi radius. Thus, the length L is size-matched 
to = {8'k/15)R^p JH m. Taking the peak density 
Po — 7.6 X 10^^ m~^, the trap frequency uj /2Tr = 52 
Hz, and molecule-molecule scattering length a,„„ ~ 1.47 
pm in the strongly interacting regime near a magnetic Fesh- 
bach resonance, we obtain Rtf — \/ ^''^fi-^O'mmPa I ttv^'-lP- ~ 
33.1 pva, L ~ 39.3 pva, and a total number of molecules 
A^o(O) ~ 3 X 10^. Assuming a 3D atom-molecule coupling 
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Eq. (IA5l l. By energy and momentum conservation, the to- 
tal dissociation energy 2fi.| A| is converted into the kinetic en- 
ergy, 2h\A\ h?ky{2mi) + h?k1/{2mi), of atom pairs 
with opposite spins and momenta. The atomic momenta 
are distributed in a narrow interval around fco = |k|(x)| = 



FIG. 2: Dissociation of a molecular condensate into distinguishable 
bosons. Apart from bosonic statistics, all parameters are the same as 
for Fig.[T] 



XzD = 4.87 X 10 ^ m^/'^s we obtain a time scale of to = 
1/[ka/]Vo(0)] - 0.918ms, where K Xzu/L^'"^ ~ 1.98 s^^ 
and giving a length scale of do — 0.854 /im, I = L/da :^ 46, 
and hence A^o(O)//^ ^ 3.1. 

Most of the numerical examples given here are for two val- 
ues of the dimensionless detuning, \S\ — 16 and \S\ = 3174, 
where 6 — Ato. Using the above timescale of ^o — 0.918 
ms, these detunings correspond to |A|/27r ~ 3.67 kHz and 
|A|/27r ~ 0.55 MHz. The second case is chosen to coincide 
with the average rf detuning of I'rf = 11 MHz used in the 
experiment of Ref. |5], where we note that u,./ corresponds to 
the absolute detuning of |A|/7r in our notation. The duration 
of dissociation of 0.33 ms used in Ref. corresponds to di- 
mensionless time of T ~ 0.36 at which our 3D results give 
Ni{t)/NQ{0) ~ 0.39 (39% conversion), in good agreement 
with ~ 43% conversion obtained in Ref. Jst]. We also note 
that the dimensionless time Tmax = corresponding 
to Eq. (|9]l is equal to r^ax — 0.41 in this example, which is 
larger than the dissociation duration of r ~ 0.36 as required 
for the applicability of the uniform treatment. 

More detailed quantitative comparison with the results of 
Ref. fs^l, especially those for the momentum distribution and 
atom-atom correlations, is complicated by the fact that the de- 
tuning Vrf was swept over 600 kHz during the experiment in 
order to spread the detected signal over many camera pixels. 
We do not model the rf sweep in the present work. Additional 
discrepancies are expected to arise due to the fact that the trap 
geometry in the experiment was not spherically symmetric as 
assumed here. 

The results of Figs. [T] and |2] are easily understood using 
the dependence of the density of atomic states on energy. 



^y2ml\A\/h. In 3D, the resonant momenta form a spheri- 
cal shell with a radius fco; in 2D they are distributed within 
a circular shell at the same radius, while in ID the resonant 
momenta are around ±fco. Since the density of atomic states 
at the dissociation energy is the largest in 3D, the number of 
dissociated atoms is larger in 3D than in 2D and ID. 

Next, we address the question of how the number of dis- 
sociated atoms changes with the dissociation energy 2/i|A|. 
Since the density of states increases with energy in 3D, a 
larger absolute detuning results in a larger number of atoms 
produced [compare the 3D curves in Figs. [Ha) and (b); also in 
Figs. 121 a) and (b)]. In ID the situation is reversed; the density 
of states decreases with energy and hence the number of atoms 
at a given time is smaller at larger absolute detuning |A|. In 
2D, the density of states is independent of the energy and we 
see no variation of the total atom number with the detuning. 
Indeed, the 2D curves in Figs. [1] and |2] are indistinguishable 
within numerical accuracy, even though the equations were 
solved with different detunings. 

As expected, the dependencies on the absolute detuning in 
ID, 2D and 3D found numerically are in agreement with the 
explicit analytic results for the atom production rate found 
from Fermi's golden rule, Eq. ( fT2] l. The atom numbers from 
the simple rate equation ( fT2] l are shown by dotted lines in Figs. 
[T] and 121 where they appear as the tangents to the respective 
numerical results at r ^ 0. 

Comparing the results of the pairing mean-field theory with 
those based on Fermi's golden rule calculations of the molec- 
ular decay rate, we see that in the initial spontaneous stage of 
dissociation the bosonic and fermionic results are very simi- 
lar Past the spontaneous regime, the dynamics of dissociation 
becomes affected by either Pauli blocking or bosonic stimu- 
lation, depending on the statistics of the dissociated atoms. 
Accordingly, the fermionic results for the total atom number 
saturate faster and lie below the curve corresponding to the 
simple rate equation, while the bosonic results exhibit expo- 
nential growth due to Bose enhancement and remain above 
the rate equation curve. 

In the bosonic case, the ultimate limit to the growth in atom 
number is set by the total (finite) initial number of molecules 
as the condensate can be entirely depleted. For fermions, how- 
ever, the Pauli exclusion principle - depending on dimension- 
ality of space and absolute detuning - can take effect before 
there is any significant depletion of the molecular condensate. 
In this sense the molecular depletion is a more important fac- 
tor in the bosonic case than in the fermionic case. This sit- 
uation applies to all curves in Fig. [Ha) and to ID and 2D 
results in Fig. [Tfb). In all these fermionic cases the conver- 
sion efficiency is less than 8% during the entire simulation 
time and the molecular depletion is negligible. Accordingly, 
the dynamics of dissociation follows closely the predictions 
obtained within the undepleted molecular field approximation 
of Ref. |50]. This is shown in Fig. [3ja), where we present 
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FIG. 3: Same as in Fig. [TJa) and (b) in thie 3D case, except that the re- 
sults of the pairing mean-field theory (solid lines) are compared with 
those obtained within the undepleted molecular field approximation 
(dashed line). 



2. Parametrization at large dissociation energy 



Once the equations of motion are written in dimensionless 
form it is clear the system can be described via just two param- 
eters — the dimensionless detuning 5 and the initial number 
of molecules Nq ~ 0^. The quantization volume drops out of 
the final results if we are interested in fractional populations 
or normalized densities rather than their absolute values. 

For large dissociation energy (large absolute detuning, 
I (5 1 ^ 1), the parametrization of the system can be further 
simplified and reduced to just one parameter which is a func- 
tion of iVo and \5\ ^231]. Here, we show this for the case of 
distinguishable atoms and note that the same arguments apply 
to the indistinguishable case. As shown in Appendix ICl the 
parameter in question originates from the approximate form 
of Eqs. ( |23] ) in the continuous limit, which gives the follow- 
ing equation for the fractional amplitude 



df{r) 
dr 



dSkm{Sk,T). 



(26) 



Here, the dimensionless parameter T is defined via 



T 



2No 



h^K^[D{h\A\)]^' 



(27) 



the comparison in the 3D case, where the discrepancy is the 
largest. Owing to the different dependencies of the density of 
states on the absolute detuning in different dimensions, the un- 
depleted molecular field approximation works better at small 
absolute detunings in 3D and at large detunings in ID. For the 
same reason, the molecular depletion generally is less impor- 
tant in ID than in 3D. 

If, on the other hand, the number of available atomic states 
in the spherical shell around ko [see Eq. ( l35T l for the width of 
the shell] is comparable to or larger than the total initial num- 
ber of molecules, then the dynamics of dissociation is domi- 
nated by the molecular field depletion. In this case, the num- 
ber of molecules can decrease substantially before the popula- 
tion of individual atomic modes experiences any Pauli block- 
ing. This is a typical situation in 3D at very large absolute 
detunings \6\ when there is a large number of states available 
for occupancy. Even if the average occupation of each of these 
states is smaller than one, the total average number of atoms 
can be quite large and constitute a large fraction of the initial 
number of molecules. Accordingly, the predictions of the un- 
depleted molecular field approximation remain valid only for 
very short time as seen in Fig. |3jb). In the longer time limit, 
the undepleted molecular field approximation leads to an un- 
physical result that the fractional atomic population becomes 
larger than one, as seen in Fig. [3lb). This situation is simi- 
lar to the well-known behavior in the bosonic case (see Fig.|5] 
below). 



and coincides with the parameter F introduced in Ref. lEsIl 
(where we note that N — 2A'o); this is the only parameter 
that characterizes the system at large detunings. Using Eq. 
dASI l, T can be written down explicitly for ID, 2D, and 3D as 
follows: 

r Sir^m^/P, (ID), 
T = J 327^2iV2/^^ (2D), (28) 
[ 327r^Ni/im), (3D). 

This approximation breaks down for small absolute detun- 
ings where the dissociation predominantly populates a range 
of atomic momenta close to zero and the density of states here 
can no longer be approximated as flat. As a result the system 
is parameterized by two variables T and \5\, or equivalently 
in terms of the original pair Nq and \S\, without the need to 
introduce T. 

In Fig. m we show the comparison between the results ob- 
tained with and without the use of the flat density of states 
approximation. In this example, the absolute dimensionless 
detuning is \5\ = 1 and we see significant discrepancy be- 
tween the two curves. The discrepancy increases further for 
smaller detunings, while it becomes negligible for detunings 
l-^l > 16. 
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FIG. 4: Comparison of the results obtained withiout (solid line) and 
with (dash-dotted line) the flat density-of-states approximation in 3D. 
In both cases we set |^| = 1 and No/l'^ = 3.1, corresponding to 
T = 3 X 10*. For the dash-dotted line, only the value of T matters, 
and the same curve can be obtained with different individual values 
of \5\ and No/l'^ as long as they result in the same T, Eq. i2Si . 



FIG. 5: Fractional atomic population Ni{t) /2No{0) in 3D ob- 
tained using the present pairing mean-field theory (PMFT) and first- 
principles (exact) simulations of the same system using the positive- 
P representation method (+P) i24ll . The dashed line corresponds to 
the undepleted molecular field approximation. The parameter values 
are: \S\ = 4 and iVo(O)//^ = 0.988. 



B. Indistinguishable atoms 

In the case of dissociation into indistinguishable atoms our 
treatment corresponds to the Hamiltonian 

-i^ ^^{blai^ai-k - _k«i,k^o)' (29) 

where again n = xd/L^^"^ and xn {D — 1, 2, 3) are given 
by Eqs. (Q, CJ, and dill 

The analysis of this system within the PMFT is essentially 
the same as in the distinguishable case, Eq. ( l23T l. The only 
difference is in the equation for /(r), which now reads as 

while the normal and anomalous populations are 

nk(t) = {niMit)) = (alkW«i,k(i)), (31) 
mk(t) = (mk(t)> = (ai,k(i)ai,-kW>. (32) 

In Fig. |5]we show a comparison between the present pairing 
mean-field theory and first-principles (exact) simulations of 
the same system using the positive-P representation method 
ll24i] . As we see, the pairing mean-field theory compares well 
with the exact results for dissociation durations correspond- 
ing to less than 50% conversion at which stage the discrep- 
ancy reaches ~ 20%. This is much better than the results 
obtained using the undepleted molecular field approximation 
lH^l, shown by the dashed line. 

We note, however, that the comparison between the PMFT 
results and those of the positive-P method are not entirely 
equivalent as the positive-P method takes into account the 
possibility of a dynamical population of the initially unoccu- 
pied molecular modes with non-zero momenta. As has been 
shown in Ref. |24], this process becomes a sizeable effect 
as time progresses and occurs due to "rogue" association in 



which a small fraction of newly formed atoms convert back 
into molecules with nonzero momenta. The present version 
of the pairing mean-field scheme only treats the molecular 
condensate mode foo [see Eq. ( fTST ll without allowing for pop- 
ulation of the non-condensate modes 6k (k ^ 0) present 
in the original Hamiltonian ( flTl ). The performance of the 
PMFT can be improved (at the cost of increased computa- 
tional complexity) by incorporating the mean-field dynamics 
of all non-condensate modes {h\^{t)) Pk{t), irrespective 
of their initial population. In this manner one can extend 
the present treatment from uniform to physically more realis- 
tic nonuniform systems corresponding to spatially inhomoge- 
neous trapped condensates. This will be considered in future 
work. 



V. MOMENTUM DISTRIBUTION AND ATOM-ATOM 
CORRELATIONS 

A. Atomic momentum distribution 

In this section we analyze the momentum distribution of 
the dissociated atoms, nk(t), which is the same for both spin 
states. In Fig. |6ja) we plot the momentum distribution ver- 
sus time for the case of fermionic atom statistics. Due to the 
spherical symmetry in 3D, we only show the radial depen- 
dence on the absolute momentum k = |k|, scaled with respect 
to the resonant momentum ko = ■\/2mi| A|/fi,. 

In Fig.l^b) we show examples of snapshots of the momen- 
tum distribution at times r = 1 and t — 3. The respective ID 
and 2D results are very similar to the 3D results and are not 
shown. The maximum deviation from the 3D result is in the 
ID case and is less than 2% at r = 3. 

In Fig.|6jc) we plot the occupation of the resonant mode fco 
as a function of time. The dashed line is the analytic result 
obtained using the undepleted molecular field approximation 
is given by sin'^(T), and is applicable in all dimensions jSOil . 
Again the ID and 2D results of the present pairing mean-field 
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FIG. 6: Atomic momentum distribution in each spin state for disso- 
ciation into fermionic atoms for a dimensionless absolute detuning 
\5\ = \A\to = 16 and No{0)/l° = 3.1. Only the range of momenta 
whose population is nonvanishing is shown, (a) nk(T) as a function 
of scaled time r = t/to and scaled absolute momentum k/ko in 3D. 
(b) Time slices of the momentum distribution at r = 1 and t — 3. (c) 
Temporal population of the resonant momentum mode with k — ko. 
The dashed line represents the respective analytic solution given by 
sin'^(r) in the undepleted molecular field approximation L5Q1 . 



theory are omitted for clarity — in the ID case the result is 
almost indistinguishable from the dashed line, while in 2D it 
is intermediate between the dashed line and the 3D solid line. 
The larger discrepancy of the 3D result from the approximate 
analytic solution is explained by the fact that molecular de- 
pletion is a more significant effect in 3D. This is better seen in 
Fig.|7]for the significantly larger absolute detuning \d\, where 
the molecular depletion is more important and hence the dis- 
crepancy between the pairing mean-field result and that of the 
undepleted molecular field approximation is larger. 

Figures [8] and |9] are the results for the same parameters as 
in Figs.|6]and|2]respectively except that they are for bosonic 
atom statistics. Compared to the fermionic case of Fig.|6l the 
dissociation dynamics in Figs. [8] are not affected by the Pauli 
exclusion principle and we see larger discrepancies between 
the ID, 2D and 3D results. The disagreement between the un- 
depleted molecular field approximation and that of the present 
PMFT is again largest in the 3D case due to the larger density 
of states which results in faster depletion. 

It is instructive to use the above numerical results to de- 
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FIG. 7: Same as in Fig. IHexcept for \5\ = 3174. In (c) we show the 
temporal population of the resonant mode fco in ID, 2D, and 3D. The 
analytic result of sin^(r) corresponding to the undepleted molecular 
field approximation is practically identical to the ID curve shown. 



termine the width of the shell of the dissociation sphere and 
compare it with the predictions of much simpler analytic ap- 
proaches. By doing this we can get a better understanding 
of the validity of these approximate methods. For exam- 
ple, the simplest estimate of the width of the dissociation 
sphere can be obtained from energy-time uncertainty consid- 
erations. More specifically, for short duration of dissociation 
t (Ft ^ 1) during which the dissociation is in the sponta- 
neous regime and is independent of atom statistics, the width 
of the dissociation sphere can be obtained through the rela- 
tion AEt ^ h. Here, the energy uncertainty is evaluated in 
terms of the momentum uncertainty Afc around the resonant 



momentum feg — \/2mi\A\/h: 



AE = 



2mi 



2mi 



h^koAk 
mi ' 



(33) 



where we assumed Afc/fco <C 1. According to this, the mo- 
mentum uncertainty Ak which gives the width of the dissoci- 
ation sphere is given by 



Ak 

ko 



1 



1 



hk^t 



2\A\t 2\6\t' 



(34) 



which is inversely proportional to the duration of dissocia- 
tion t as expected and also to kg (or |A|), implying that for 




FIG. 8: Same as in Fig. [6](|(5| = 16) except for bosonic statistics of 
the atoms. The time slices of the momentum distribution in (b) are 
for r = 1 and r = 2. In (c), we show the temporal population of the 
resonant mode ko in ID, 2D, and 3D (solid line), together with the 
respective analytic result, sinh'^ (r) (dashed line), obtained using the 
undepleted molecular field approximation (50ll . 



larger detuning | A| the width of the sphere is narrower. These 
conclusions are clearly seen from our numerical results for 
T < 1, shown in Figs. |6]-|9] For example, at t ~ 1 we 
obtain [from Eq. dSli] Afc/fco ~ 0.03 for \S\ = 16 and 
A/c/fco ~ 1.6 X 10""* for \S\ = 3174, in good agreement 
with the numerical results of Figs.|6jb)-[9jb). 

For longer duration of dissociation, when the dynamics is 
strongly affected by either Bose stimulation or Pauli block- 
ing, the width of the dissociation sphere can be estimated 
(see AppendixiDli using a slightly more involved approximate 
approach based on analytic solutions within the undepleted 
molecular BEC approximation ll24ll50ll . In the fermionic case, 
where the momentum spectrum develops oscillatory peaks, 
the width in question is for the envelope function. Apart from 
this distinction between the fermionic and bosonic results, the 
width of the dissociation spherical shell is the same in the two 
cases and is given by 



Afc ^ mxD^/N^ 
ko 2hkl 



m' 



(35) 



FIG. 9: Same as in Fig. |7](|(5| = 3174) except for bosonic atom 
statistics. 



distribution in this regime (t > 2 in most of our examples) 
settles into a shape with an approximately constant (envelope) 
width and residual oscillations [see Figs. |6]-|9)- For \5\ = 16 
and \6\ — 3174, the above expression gives, respectively, 
Afc/fco = 0.016 and Afc/fco = 0.79 x lO""* which are again 
in good agreement with the numerical results. 



B. Atom-atom correlations 

1. Distinguishable case 

In this section we study the correlation between the occu- 
pation number fluctuations for atoms in opposite spin states 



Gi2(k,k',t) = 



(Ani,k(0An2,k'W) 
(ni,kW)i/2(n2,k'W>'/' 



(36) 



where Ahjy^{t) — nj\^{t) 



,k(0) is the fluctuation. 



The definition of 6*12 is different from the more conventional 
Glauber second-order correlation function 06311 



This is still inversely proportional to the dimensionless de- 
tuning |(5|, but no longer depends on time as the momentum 
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in that G12 is normalized to the square root of the product of 
the occupation numbers (rather than to the direct product) and 
the contribution from uncorrected statistics is subtracted from 
it so that absence of correlation corresponds to G12 = 0. The 
above definition of G12 has been used in Ref. lUt] for quan- 
tifying the strength of correlation between the atoms in two 
hyperfine states of ^"K atoms created through molecular dis- 
sociation, and it is particularly convenient for fermionic atom 
statistics since in this case 6*12 is bounded between and 1. 

(2) 

The correlation functions G12 and are related by a simple 
relationship. 



5r2^(k,k',t) = l 



Gl2(k,k',i) 



(ni,k(i)>^/'("2,k'W>^/' 



(38) 



Using Wick's theorem to factore higher-order correlation 
functions and express them in terms of the normal and anoma- 
lous populations, (aj j.(t)aj,k' (^)) and (ai^k(i)%,k' (i)), we 
find, in particular, that the correlation between the atoms with 
equal but opposite momenta, k' — — k, is given by 



Gl2(k,-k,t) 



nk(t) 



= l±71kW. 



(39) 



Here, the upper (lower) sign refers to bosonic (fermionic) 
statistics, and we have used the fact that in the present pairing- 
mean field theory any second-order moments other than those 
in Eqs. ( l20b and ( |2TI ) are zero. In addition, we have used 
Eq. ( l25T l for expressing the anomalous occupation m\^{t) via 
the normal occupation number n\^{t). This gives the maxi- 
mum possible value of the anomalous moment rrik and hence 
the maximum strength of the second-order correlation func- 
tion between the atoms with opposite momenta ll 24l Isoll . 
Therefore the results of the present (approximate) PMFT for 
correlation functions give the respective upper bounds to their 
exact values. 

For any k' ^ — k, on the other hand, we have 



Gi2(k,k',t) = 0, k'^-k 



(40) 



which corresponds to absence of any correlation. 

For completeness, we also give the results for the correla- 
tion between the occupation fluctuations for the atoms in the 
same spin state: 

For the present pairing mean-field theory, this is given by 
Gn(k,k',i) - G22(k,k',t) = I ^""^J'^^ ^' (42) 

where again the upper (lower) sign refers to bosonic 
(fermionic) statistics. Here, the case with k' 7^ k implies un- 
correlated statistics, which in terms of Glauber's second-order 
correlation function corresponds to 5^ j^'^ (k, k' 7^ k, i) = 1. 
Similarly, Gii(k, k, t) = JT-kl^) in the bosonic case corre- 
sponds to the level of correlations between thermally bunched 



(2) 

atoms, gj^j^ (k, k, i) = 2. Finally, Gii(k, k, i) = —ni^{t) in 
the fermionic case corresponds to g\j' (k, k, =0 implying 
perfect antibunching, which is a simple consequence of the 
Pauli exclusion principle. 

The result of Eq. ( [39] l for the opposite spin atoms corre- 
sponds to 



5if(k,-k,t) = l 



l±»kft) 
nk(i) 



(43) 



Thus, strong correlation between the atoms with equal but op- 
posite momenta corresponds to (k, —k, i) — l/nk(0 > 1 
in the fermionic case, with rtk(i) < 1 due to Pauli blocking, 

(2) 

and to (k, — k, t) = 2 + l/nk(i) > 2 in the bosonic case, 
corresponding to super-thermal bunching. 

Since Wick's factorization scheme as employed here is 
equally valid for the pairing mean-field theory and for the 
treatment of Ref. |50] with the undepleted molecular field ap- 
proximation, all results for second-order correlation functions 
obtained in Ref. fsoll remain valid here. The only difference is 
that the analytic solutions for nk(i) fsdl are now replaced by 
the numerical solution of Eqs. (|23]) . This takes into account 
molecular depletion and therefore the present results are ex- 
pected to be more accurate. 



2. Indistinguishable case 

In the case of molecular dissociation into bosonic atom 
pairs in the same spin state, the considerations and analysis 
of the previous subsection remain valid, except that the role 
of Gi2(k, — k, t) is now taken by Gii(k. — k, t). We note, 
however, that for indistinguishable atoms the conveniences 
associated with using Gii(k, k',i) - defined analogously to 
Eq. (l36T l - for characterizing the strength of correlation be- 
tween atom pairs are no longer valid. For example, in the 
present pairing mean-field theory we find that the correlation 
functions Gii(k, — k, i) and Gii(k, k, i) are both equal to 
1 + rL]f^{t), which does not distinguish between the strength 
of correlation for atom pairs with opposite momenta and the 
strength of autocorrelation signal. For this reason, we use 
Glauber's second-order correlation function 



5ii'(k,k',t) = 



(«l,kWQl,k'W°l.k' ft)Ql,k(0) 

(ni,k(t)) (ni.k'W) 



(44) 



Following the same arguments as in Ref. 1,24.1 . we obtain 
the following results 



3fi^(k,k',i) 



1, k'7^k,k'7^ 

2, k' = k 7^ 0, 

2 + l/r^kW, k' = -k, 

3 + l/nk(t), k' = k = 0, 



(45) 



which again formally coincide with those obtained using the 
undepleted molecular field approximation fl^, except that 
the mode occupation numbers rik(t) are now given by the 
numerical solutions of the present pairing mean-field theory, 
Eqs. (US and 
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results for equal but opposite momenta, p' = — 



Gi2(p,-P,<) = l±^^ 



FIG. 10: (a) Momentum-space column density np{T)/l of each spin 
component at r — 0.36 for bosonic and fermionic cases in 3D for 
\S\ = 3174 and iVo(0)/;-° = 3.1. The horizontal axis is the absolute 
value of the 2D reduced momentum p — jpj, where p ={kx,ky). 
(b) Correlation function between momentum column density fluc- 
tuations for the atoms with opposite spins and conjugate momenta, 
Gi2(p, — p,r), as a function of p = |p|, for the same parameter 
values as in (a). 



C. Correlations in column densities in 3D systems 

The correlation functions derived and discussed in the pre- 
vious section are difficult to measure experimentally. Rather 
than density functions such as n(x) being measured in the 
laboratory, column densities such as J dzn{x) are obtained 
using the technique of absorption imaging. Indeed, the cor- 
relation measurements performed in the experiments at JILA 
on dissociation of a molecular condensate of *°K2 dimers into 
fermionic atom pairs |^ were made using absorption images 
after a time-of-flight expansion. We now analyze the momen- 
tum space analog of this procedure and calculate the correla- 
tion between momentum column density fluctuations. 

First, we define the atom number operator corresponding to 
a z-integrated momentum column density, nj p = fij ^^, 
where p ={kx, ky) is the reduced 2D momentum. The aver- 
age column density rip{t) — {nj^p{t)) is found via rTp(t) = 
X]fc "-kC^) and is the same for the two spin states, j = 1, 2. 
The correlation function between the momentum column den- 
sity fluctuations in the two spin states, which we denote via 
Gi2(k, k', t), is defined as in Eq. (|36] | except that the opera- 
tors rij k are replaced by nj.p: 

{mAt)V/Hn,,p,{t)y/^ 

We emphasize that the summation over the -component has 
to be performed before taking the ensemble average. 

Using Wick's theorem as before, we obtain the following 



.it) 



(47) 



Here the plus sign stands for bosons, while the minus sign 
is for fermions, and the distinction between the bosonic 
and fermionic results is highlighted through the fact that 
Gi2{p,—p,t) > 1 for bosons and Gi2(p, —p, t) < 1 for 
fermions. For any other pair of momenta (p' ^ p), the 
correlation function is zero, implying the absence of any cor- 
relation. 

In Fig. [To] we show examples of momentum column den- 
sities np{t) and correlations between column density fluc- 
tuations Gi2(p, — p,<) for bosonic and fermionic atoms for 
\S\ ~ 3174 at r = 0.36. The dissociation at this stage is still in 
the spontaneous regime, so the difference between the column 
densities for bosonic and fermionic cases is hardly noticeable 
in Fig. [Tol a). On the other hand, the difference due to quantum 
statistics becomes much more pronounced if one examines the 
pair correlation function Gi2(p, — p, t), as shown Fig. [Tol b). 
with Gi2(p, — p, t) > 1 for bosons and Gi2(p, —p, i) < 1 
for fermions. 

The example of Fig. [TO] in the fermionic case, with 
Gi2(p, — p, t) ~ 0.9, corresponds to the experimental param- 
eters of Ref . f3\ . We note that the experimentally measured 
spatial counterpart of the same correlation function after time- 
of-flight expansion varied between 0.1 and 0.3 depending on 
the bin size used in the analysis of the imaged data. The de- 
pendence of the strength of the correlation signal on bin size 
was analyzed theoretically in Ref. |31]; However, ignoring 
variations due to binning, the experimental result was sig- 
nificantly less than our theoretically calculated value in mo- 
mentum space. This highlights an important limitation of the 
present uniform theory, which is that it only gives the upper 
bound of correlation functions |50] and is not adequate for 
correct quantitative description of realistic nonuniform sys- 
tems. This conclusion, however, does not necessarily apply 
to other quantities of interest, such as the momentum distri- 
bution or the total atom number. In these cases, the quantita- 
tive predictions of the uniform theory can adequately describe 
nonuniform systems if the uniform system is properly size- 
matched to the nonuniform one as done in Sec. IIV A 1 1 (see 
also Ref. iH). 

The reason that the PMFT gives unsatisfactory results for 
pair correlations is the mode-mixing which couples the atomic 
momenta not only in a pair-wise fashion, (k,— k), but also 
couples within the neighborhood of the partners. As a result 
the correlations are broadened, while their strength can be dra- 
matically reduced, especially in systems with strong inhomo- 
geneity. The problem of mode-mixing was analyzed in detail 
in the bosonic case in Ref. |24] using the comparison between 
uniform solutions and exact quantum simulations of nonuni- 
form systems using the stochastic positive P-representation 
method. Performing similar first-principle quantum dynami- 
cal calculations for nonuniform, muItimode/erm/o«/c systems 
is a challenging problem yet to be solved. 
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VI. CONCLUSIONS 

We have developed a simple pairing mean-field theory for 
the quantum dynamics of the dissociation of a Bose-Einstein 
condensate of diatomic molecules into its constituent parti- 
cles, covering the cases of indistinguishable bosons and dis- 
tinguishable fermions or bosons. The pairing fields for the 
atoms are introduced at the level of normal and anomalous 
moments for the atomic creation and annihilation operators; 
such an approach is necessary as the dissociation is initially 
a spontaneous process and cannot be described by ordinary 
mean field theory. We have solved the resulting equations of 
motion for all cases, and identified the role of dimensionality 
and quantum statistics on the conversion rate, and compared 
them to the simpler analytic solutions available within the un- 
depleted molecular field approximation. 

We find that molecular depletion is a more important ef- 
fect in higher dimensions and in the case of dissociation into 
bosonic atoms. Accordingly, we argue that the undepleted 
molecular approximation is more reliable for describing dis- 
sociation into fermionic atoms and systems of reduced dimen- 
sionality. This conclusion can be useful in extensions of the 
theory of dissociation to spatially inhomogeneous systems. In 
this case, the simplest treatment would be to employ the un- 
depleted molecular field approximation, while still being able 
to address the problem of mode-mixing, which is the most 
important one when seeking quantitative description of atom- 
atom correlations. 

We have also compared the results in the bosonic, indis- 
tinguishable case to the exact dynamics found by positive-P 
simulation and identified the range of validity of the pairing 
mean-field theory. In addition, we have compared our formal- 
ism to that of Jack and Pu 1,23,1. and identified the key dif- 
ferences. Finally, we have derived and calculated results for 
the non-trivial correlations that are present in the dissociated 
atoms, and compared our results to those of the experiment 
performed at JILA Jl]. 

The main limitation of this work (see also Sec. IllCb is the 
assumption of a uniform molecular condensate. This limits 
the time for which the dynamics are valid and also does not 
allow for the association of atoms into previously unoccupied 
molecular modes. The PMFT can be used to describe realistic 
nonuniform condensates provided the results are applied to a 
size-matched system. The total atom number and atomic den- 
sity distribution obtained in this way can be quantitatively rea- 
sonable, however, the same is not true if the calculated quanti- 
ties involve density-density correlation functions. In this case 
the PMFT provides upper bounds for the correlation func- 
tions. We are currently extending the PMFT to include the 
effects of spatial inhomogeneity to provide a quantitatively 
better model. 
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APPENDIX A: DECAY RATE IN THE DISTINGUISHABLE 
CASE 

The transition matrix element Vma in Eq. ( fTOl i is 

Vma = {f\H,nt\i). (Al) 

where Hint is the atom-molecule coupling term in the Hamil- 
tonian ([T]l, while the initial |i) and final |/) states describe a 
free molecule or a pair of free atoms respectively 



1 



d^xv|/T(x)e^K-|0), 



l/> = 



^ 1 1 d^^d^y 



jK-(x+y)/2 



|0). 



(A2) 



(A3) 



Here K is the total center-of-mass momentum, L is the quanti- 
zation length in each dimension, and the states are normalized 
to one. Calculating the matrix element in Eq. ( lAll ) gives 



,y ,2 Xd 
\yma\ — ■ 



(A4) 



Next, the density of two-atom states D^^' (e) at the total dis- 
sociation energy e = 2h\A\ is the same as the standard den- 
sity of single-atom states D{E) at energy E ~ e/2 evaluated 
at half the quantization volume V = L'^ for a single molecule 
(or area A ~ in 2D, or length L in ID); equivalently it 
is the same as half the density of single-atom states evaluated 
for the same quantization volume, and is given by iBTll 



aL/{2VE), (ID), 
7r(aL)V2, (2D), (A5) 
^ n{aLfVE, (3D), 



where a — {2mi/fpy^'^/2n is a constant. Combining 
Eqs. ([Tot . (IA4l i and (IA5b . and expressing the final result in 
terms of the absolute detuning | A|, we obtain Eq. (fTTl i for the 
decay rate F. 



APPENDIX B: DECAY RATE IN THE 
INDISTINGUISHABLE CASE 

The initial \i) and final |/) states describing, a free molecule 
and a pair of free indistinguishable atoms respectively are 



K> = 
l/> = 



1 



1 



d^x*T(x)e*' 



V2LD 

x*I(x)*I(y)e'^-("+y)/2|o), 



(Bl) 



(B2) 



and are again normalized to one. Applying these to the tran- 
sition matrix element Vma = {f\Hint\i), where Hint is the 
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interaction term in the Hamiltonian (|6]l, we obtain 

W^a? = (B3) 

Accordingly, the molecular decay rate, V = 
2T:\Vma?D^'^\e)/h, where D^'^\e) is given by Eq. ( lASl l, 
takes the form of Eq. ( fTsT l. 



for a uniforrn system in the undepleted molecular BEC ap- 
proximation 1241 [soli . The solution for the momentum distri- 
bution nk;(t) in the bosonic case reads as 

Mt) = -^2^ silV {jg^-^lt^ , (Dl) 

where 



APPENDIX C: PARAMETRIZATION AT LARGE 
DISSOCIATION ENERGY 



Following the arguments of Ref. 12311 we first take the con- 
tinuum limit of Eqs. ( |23] | and convert the sums over k into 
integrals. Changing from the momentum variable to the en- 
ergy Ek — h'^k'^ /2mi and then to hAk = + h/S. where we 
recall that A < 0, we can rewrite the equation for /(t) as 



df{T) 
dr 



|A| 



dAkD{hAk-fiA)mk{T), (CI) 



where D{E) is the density of states Eq. 

At large |A| the resonance condition ~ (or Ek = 
fi| A|) leads to the population of only those atomic states that 
have absolute momenta k = |k| in a narrow range around 
ko ~ ^y2ml\A\/h. The density of states in this range of 
momenta is essentially constant and therefore D{Ek) can be 
approximated by D{h\A\) and be taken out of the integral in 
Eq. ( ICll ). At the same time, the lower limit of the integral 
(— I A|) can be replaced by — oo, which leads to the following 
result in terms of dimensionless parameters 



df{T) hD{h\A\) 



dr 



dSkm{5k,T), (C2) 



A. 



XD 

he 

2m 



hkl 

2m' 



(D2) 
(D3) 

From Eq. jDll l we see that modes with — A^ > expe- 
rience Bose enhancement and grow exponentially with time, 
and the modes with — A\ < oscillate at the spontaneous 
noise level. The absolute momenta of the exponentially grow- 
ing modes lie near the resonant momentum /cq, and therefore 
we can use the condition g^ — A^ = to define the width 
of the dissociation sphere. First we write k ~ k^ ^ 5k and 
assume for simplicity that /jq is large enough so that 5k <C fco- 
Then the above condition can be approximated by 



(D4) 



This can be solved for 5k and used to define the width Ak 
5k/ 2 of the dissociation sphere as 



Ak 

ko 



mg _ mxDy/No/L^ 



2hkl 



2hkl 



(D5) 



where 5k = toAk/h and ■m{5k,T) is a continuous function 
corresponding to mkir) after making the variable changes. 

The coefficient in front of the integral in Eq. ( IC2b motivates 
the introduction of a dimensionless parameter T defined via 

hD{h\A\)/{(3^to) = V27T and therefore 



T = 



2l5ttl 



2iVn 



h^{D(h\A\)f h^K^\D(h\A\)Y' 



(C3) 



The reason for defining it as half the full-width 5k is to make 
Afc closer in definition to the rms width around feg- 

For the case of fermionic statistics a similar analytic solu- 
tion exists l,5Q1 



nk(i) = 



r 



A?. 



sm 



(D6) 



The final explicit result for T in ID, 2D, and 3D is given in 
Eq. (|28]l. 



APPENDIX D: WIDTH OF THE DISSOCIATION SPHERE 

To estimate the width of the shell of the dissociation sphere 
beyond the spontaneous regime we use the analytic solutions 



where g and A^ are the same as in Eqs. ( ID2I ) and ( ID3I ). In 
this case the solutions for all k are oscillatory, with popula- 
tions of individual modes not exceeding one due to the Pauli 
exclusion principle. The width of the dissociation sphere can 
now be defined as the width of the envelope function given 
by the Lorentzian g"^ /[g"^ + A^]. Rewriting the width of this 
Lorentzian in terms of the momentum width Afc = 5k/2, we 
obtain the same expression as in Eq. ( ID5b . 
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